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Abstract. Competing phases or interactions in complex many-particle systems 
can result in free energy barriers that strongly suppress thermal equilibration. Here 
we discuss how extended ensemble Monte Carlo simulations can be used to study 
the equilibrium behavior of such systems. Special focus will be given to a recently 
developed adaptive Monte Carlo technique that is capable to explore and overcome 
the entropic barriers which cause the slow-down. We discuss this technique in the 
context of broad-histogram Monte Carlo algorithms as well as its application to 
replica-exchange methods such as parallel tempering. We briefly discuss a number 
of examples including low-temperature states of magnetic systems with competing 
interactions and dense liquids. 

1 Introduction 

The free energy landscapes of complex many-body systems with competing 
phases or interactions are often characterized by many local minima that are 
separated by entropic barriers. The simulation of such systems with conven- 
tional Monte Carlo [1] or molecular dynamics [2] methods is slowed down by 
long relaxation times due to the suppressed tunneling through these barriers. 
While at second order phase transitions this slow-down can be overcome by 
improved updating techniques, such as cluster updates [3,4], this is not the 
case for systems which undergo a first-order phase transition or for systems 
that exhibit frustration or disorder. For these systems one instead aims at im- 
proving the way that relatively simple, local updates are accepted or rejected 
in the sampling process by introducing artificial statistical ensembles such 
that tunneling times through barriers are reduced and autocorrelation effects 
minimized. In the following we discuss recently developed techniques to find 
statistical ensembles that optimize the performance of Monte Carlo sampling. 
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first in tlie context of broad-histogram Monte Carlo algoritliins and tlien out- 
line how these methods can be applied in the context of replica-exchange or 
parallel-tempering algorithms. 



2 Extended Ensemble Methods 



Let us consider a first-order phase transition, such as in a two-dimensional 
g-state Potts model [5] with a Hamiltonian 

H = -JY,6.^.^, (1) 

(i,0) 

where the spins cji can take the integer values 1, . . . ^q. For q> A this model 
exhibits a first-order phasc^ transition, accompanied by exponential slowing 
down of conventional local-update algorithms. The exponential slow-down 
is caused by the free-energy barrier between the two coexisting meta-stable 
states at the first-order phase transition. 

This barrier can be quantified by considering the energy histogram 

ffca„o„ical(£^) CX .g(S)PBoltzma„n (^) = g{E) exp(-/3S) , (2) 

which is the probability of encountering a configuration with energy E during 
the Monte Carlo simulation. The density of states is given by 

g{E) = Y,5EMc), (3) 

c 

where the sum runs over all configurations c. Away from first-order phase 
transitions, ifcanonicai(-E') has approximately Gaussian shape, centered around 
the mean energy. At first-order phase transitions, where the energy jumps dis- 
contimiously, the histogram i?canonicai(^') develops a double-peak structure. 
The minimum of ifcanonicaK-E) between these two peaks, which the simulation 
has to cross in order to go from one phase to the other, becomes exponen- 
tially small upon increasing the system size. This leads to exponentially large 
tunneling and autocorrelation times. 

This tunneling problem at first-order phase transitions can be aleviated 
by extended ensemble techniques which aim at broadening the sampled en- 
ergy space. Instead of weighting a configuration c with energy E = E{c) us- 
ing the Boltzmann weight Psoitzmann (-2^) = exp{—f3E) more general weights 

-^extended 

(E) are introduced which define the extended ensemble. The config- 
uration space is explored by generating a Markov chain of configurations 

Ci ^ C2 ^ ■ . . ^ Ci ^ d+i ^ ... , (4) 

where a move from configuration Ci to C2 is accepted with probability 

Wcxtcndcd(-E2) 



-Pacc(ci ^ C2) = min 



1 ^ 

'^(ci) 



' Wextended(-El) 



(5) 
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In general, the extended weights are defined in a single coordinate, such as 
the energy, thereby projecting the random walk in configuration space to a 
random walk in energy space 

Ei^E{ci)^E2^ ...^Ei^Ei+i^ ... . (6) 

For this random walk in energy space a histogram can be recorded which has 
the characteristic form 

-f^oxtcndcd(-E) OC .g(-E)VFoxtondod(^^) , (7) 

where the density of states g{E) is fixed for the simulated system. 

One choice of generalized weights is the multicanonical ensemble [6] where 

the weight of a configuration c is defined as Wmuiticanonicai(c) oc l/g{E{c)). 
The multicanonical ensemble then leads to a flat histogram in energy space 

^fmulticanonical(i^) OC 5f(i^) Wmulticanonical(-E') = 9{E)^^ = COnst. (8) 

removing the exponentially small minimum in the canonical distribution. Af- 
ter performing a simulation, measurements in the multicanonical ensemble 
are reweighted by a factor WBoUzmann(-E)/VFinuiticanonicai(-E) to obtain aver- 
ages in the canonical ensemble. 

Since the density of states and thus the multicanonical weights are not 
known initially, a scalable algorithm to estimate these quantities is needed. 
The Wang-Landau algorithm [7] is a simple but efficient iterative method 
to obtain good approximates of the density of states g{E) and the multi- 
canonical weights Wmuiticanonicai(-E') c< l/g{E). Besidcs overcoming the expo- 
nentially suppressed tunneling problem at first-order phase transitions, the 
Wang-Landau algorithm calculates the generalized density of states g{E) in 
an iterative procedure. The knowledge of the density of states g{E) then 
allows the direct calculation of the free energy from the partition function, 
Z = YIe 9(-E)'^~^^ ■ ^^(^ internal energy, entropy, specific heat and other 
thermal properties are easily obtained as well, by differentiating the free en- 
ergy. By additionally measuring the averages A{E) of other observables A as 
a function of the energy E, thermal expectation values can be obtained at 
arbitrary inverse temperatures (3 by performing just a single simulation; 

(Am = ^EA{E)g{E)e-^- 



3 Mcirkov chains and random walks in energy space 

The multicanonical ensemble and Wang-Landau algorithm both project a 
random walk in high-dimensional configuration space to a onc^-dimcnsional 
random walk in energy space where all energy levels are sampled equally 
often. It is important to note that the random walk in configuration space. 
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Fig. 1. Local diffusivity D(E) of a random walk sampling a flat histogram in energy 
space for the two-dimensional ferromagnetic Ising model with TV = 20 x 20 spins. 
The local diffusivity strongly depends on the energy with a strong suppression 
around the critical energy Ec ~ —1.41 and the ground-state energy Eo — ~2N. 

equation (4), is a biased Markovian random walk, while the projected random 
walk in energy space, equation (6), is non-Markovian, as memory is stored 
in the configuration. This becomes evident as the system approaches a phase 
transition in the random walk: While the energy no longer reflects from which 
side the phase transition is approached, the current configuration may still 
reflect the actual phase the system has visited most recently. In the case of 
the ferromagnetic Ising model, the order parameter for a given configuration 
at the critical energy Ec ^ —IAIN (in two space dimensions) will reveal 
whether the system is approaching the transition from the magnetically or- 
dered (lower energies) or disordered side (higher energies). 

This loss of information in the projection of the random walk in configura- 
tion space has important consequences for the random walk in energy space. 
Most strikingly, the local diffusivity of a random walker in energy space, 
which for a diffusion time to is defined as 



is not independent of the location in energy space. This is illustrated in Fig. 1 
for the two-dimensional Ising ferromagnet. Below the phase transition around 
Ec ~ —1.417V a clear minimum evolves in the local diffusivity. In this region 
large ordered domains are formed and by moving the domain boundaries 
through local spin flips only small energy changes are induced resulting in a 
suppressed local diffusivity in energy space. 



D{E, to) = {{E{t) - E{t + tD)f)/tD 



(10) 
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Because of the strong energy dependence of the local diffusivity the sim- 
ulation of a niulticanonical ensemble sampling all energy levels equally often 
turns out to be suboptimal [8] . The performance of flat-histogram algorithms 
can be quantified for classical spin models such as the ferromagnet where the 
number of energy levels is given by [— 2A^, +2A'^] and thereby scales with the 
number of spins N in the system. When measuring the typical round-trip 
time between the two extremal energies for multicanonical simulations, these 
round-trip times r are found to scale like 



showing a power-law deviation from the A^^-scaling behavior of a completely 
unbiased random walk. Here L is the linear system size and z a critical ex- 
ponent describing the slow-down of a multicanonical simulation in the prox- 
imity of a phase transition [8,9]. The value of z strongly depends on the 
simulated model and the dimensionality of the problem. In two dimensions 
the exponent increases from z = 0.74 for the ferromagnet as one introduces 
competing interactions leading to frustration and disorder. The exponent be- 
comes = 1.73 for the two-dimensional fully frustrated Ising model which is 
defined by a Hamiltonian 



where the spins around any given plaquette of four spins are frustrated, e.g. 
by choosing the couplings along three bonds to be Jij = — 1 (ferromagnetic) 
and Jij = +1 (antiferromagnetic) for the remaining bond. For the spin glass 
, where the couplings Jij are randomly chosen to be -1-1 or -1, exponential 
scaling {z = oo) is found [8,10]. Increasing the spatial dimension of the fer- 
romagnet the exponent is found to decrease as z ^ 1.81,0.74 and 0.44 for 
dimension d = 1,2 and 3 and z vanishes for the mean- field model in the limit 
of infinite dimensions [9]. 

4 Optimized ensembles 

The observed polynomial slow-down of the multicanonical ensemble poses 
the question whether for a given model there is an optimal choice of the his- 
togram i?optimai(^^) and corresponding weights lVoptimai(^^), which ehminates 
the slow-down. To address this question an adaptive feedback algorithm has 
recently been introduced that iteratively improves the weights in an extended 
ensemble simulation leading to further improvements in the efficiency of the 
algorithm by several orders of magnitude [11]. The scaling for the optimized 
ensemble is found to scale like 0([A'' In A'']^) thereby reproducing the behavior 
of an unbiased Markovian random walk up to a logarithmic correction. 

At the heart of the algorithm lies the idea to maximize a current j of 
walkers that move from the lowest energy level, to the highest energy 



T ~ 



(11) 




(12) 
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level, E+, or vice versa, in an extended ensemble simulation by varying the 
weights T'Fcxtcndcd(£')- To measure the current a label is added to the walker 
that indicates which of the two extremal energies the walker has visited most 
recently. The two extrema act as "reflecting" and "absorbing" boundaries for 
the labeled walker: e.g., if the label is plus, a visit to does not change the 
label, so this is a "reflecting" boundary. However, a visit to E- does change 
the label, so the plus walker is absorbed at that boundary. The behavior of 
the labeled walker is not affected by its label except when it visits one of the 
extrema and the label changes. 

For the random walk in energy space, two histograms are recorded, H+{E) 
and H^{E), which for sufficiently long simulations converge to steady-state 
distributions which satisfy H+{E) + H_{E) = H{E) = W{E)g{E). For each 
energy level the fraction of random walkers which have label "plus" is then 
given by f{E) = H+{E)/H{E). The above-discussed boundary conditions 
dictate f {E.) = and f{E+) = 1. 

The steady-state current to first-order in the derivative is 

j = D{E)H{E)^, (13) 

where D{E) is the walker's diffusivity at energy E. There is no current if 
f{E) is constant, since this corresponds to the equilibrium state. Therefore 
the current is to leading order proportional to df /dE. Rearranging the above 
equation and integrating on both sides, noting that j is a constant and / 
runs from to 1, one obtains 

1 ^ dE 

j Je_ D{E)H{E) • ^ ' 

To maximize the current and thus the round-trip rate, this integral must be 
minimized. However, there is a constraint: H{E) is a probability distribu- 
tion and must remain normalized which can be enforced with a Lagrange 
multiplier: 



To minimize this integrand, the ensemble, that is the weights W{E) and 
thus the histogram H{E) are varied. At this point it is assumed that the 
dependence of D{E) on the weights can be neglected. 

The optimal histogram, iJoptimai(S), which minimizes the above integrand 
and thereby maximizes the current j is then found to be 

HoptimAE) oc ;y=py- (16) 

Thus for the optimal ensemble, the probability distribution of sampled en- 
ergy levels is simply inversely proportional to the square root of the local 
diffusivity. 
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Fig. 2. Optimized histograms for the two-dimensional ferromagnetic Ising model 
for various system sizes. After the feedback of the local diffusivity a peak evolves 
near the critical energy of the transition Ec ~ —1.41 A^. The feedback thereby shifts 
additional resources towards the bottleneck of the simulation which are identified 
by a suppressed local diffusivity. 

The optimal histogram can be approximated in a feedback loop of the 
form 

• Start with some trial weights W{E), e.g. W{E) = l/g{E). 

• Repeat 

— Reset the histograms H{E) = H+{E) = H^{E) = 0. 

— Simulate the system with the current weights for N sweeps: 

* Updates are accepted with probablity min[l, W{E')/W{E)]. 

* Record the histograms (E) and iJ_ (E) . 

— From the recorded histogram an estimate of the local diffusivity is 
obtained as 



where the fraction f{E) is given by f{E) = H+{E)/H{E) and H{E) 
is the histogram H{E) = H+{E) + H^{E). 
— Define new weights as 



D{E) oc 



1 




Increase the number of sweeps for the next iteration 



sweeps ^ -^-^^ sweeps- 
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• Stop once the histogram H{E) has converged. 

The implementation of this feedback algorithm requires to change only a 
few lines of code in the original local-update algorithm for the Ising model. 
Some additional remarks are useful: 

1. In contrast to the Wang-Landau algorithm, the weights W{E) are mod- 
ified only after a batch of A^sweeps sweeps, thereby ensuring detailed bal- 
ance between successive moves at all times. 

2. The initial value of sweeps iVgwccps should be chosen large enough that 
a couple of round trips are recorded, thereby ensuring that steady state 
data for H+{E) and H-{E) are measured. 

3. The derivative df/dE can be determined by a linear regression, where 
the number of regression points is flexible. Initial batches with the lim- 
ited statistics of only a few round trips may require a larger number of 
regression points than subsequent batches with smaller round-trip times 
and better statistics. 

4. Similar to the multicanonical ensemble, the weights W{E) can become 
very large and storing the logaritms may be advantageous. The reweight- 
ing then becomes In Woptimizcd(£^) = lnW{E) + [In ^ - lnH{E)]/2 . 

At the end of the simulation, the density of states can be estimated from 
the recorded histogram as g{E) = -ffoptimized(-E)/Woptimized(-E) and normal- 
ized as described above. 

Fig. 2 shows the optimized histogram for the two-dimensional ferromag- 
netic Ising model. The optimized histogram is no longer flat, but a peak 
evolves at the critical region aroimd Ec ~ —1.41 N of the transition. The 
feedback of the local diffusivity reallocates resources towards the bottlenecks 
of the simulation which have been identified by a suppressed local diffusivity. 

When analyzing the scaling of round-trip times for the optimized ensemble 
one finds a considerable speedup: The power-law slow-down of round-trip 
times for the flat-histogram ensemble 0{N^L^) is reduced to 0{[N In N]'^) 
for the optimized ensemble, e.g. there is only a logarithmic correction to the 
scaling of a completely unbiased random walk with 0(A^^)-scaling. For the 
two-dimensional fully frustrated Ising model the scaling of round-trip times 
is shown in Fig. 3. This scaling improvement results in a speedup by a nearly 
two orders of magnitude already for a system with some 128 x 128 spins. 

5 Simulation of dense fluids 

Extended ensembles cannot only be defined as a function of energy, but in 
arbitrary reaction coordinates R onto which a random walk in configuration 
space can be projected. The generalized weights in these reaction coordinates 



Ensemble optimization techniques 



9 



100 - 



flat-histogram 
ensemble 



100 




0.1 



0.1 



8 



16 



32 



64 



128 



1/2 



Fig. 3. Scaling of round-trip times for a random walk in energy space sampling a flat 
histogram (open squares) and the optimized histogram (solid circles) for the two- 
dimensional fully frustrated Ising model. While for the multicanonical simulation a 
power-law slow-down of the round-trip times 0{N^L^) is observed, the round-trip 
times for the optimized ensemble scale like 0([A'^lnA'^]^) thereby approaching the 
ideal 0(A'^'^)-scaling of an unbiased Markovian random walk up to a logarithmic 
correction. 

Wextcndcd(-R) are then used to bias the random walk along the reaction coor- 
dinate by accepting moves from a configuration ci with reaction coordinate 
Ri to a configuration C2 with reaction coordinate i?2 with probabihty 



The generalized weights Wextended(-R) can again be chosen in such a way that 
similar to a multicanonical simulation a flat histogram is sampled along the 
reaction coordinate by setting the weights to be inversely proportional to the 
density of states defined in the reaction coordinates, that is Wcxtcndod(-R) oc 



An optimal choice of weights can be found by measuring the local diffusiv- 
ity of a random walk along the reaction coordinates and by applying the feed- 
back method to shift weight towards the bottlenecks in the simulation. This 
generalized ensemble optimization approach has recently been illustrated for 
the simulation of dense Lennard-Jones fiuids close to the vapor-liquid equi- 
librium [12]. The interaction between particles in the fluid is described by a 
pairwise Lennard-Jones potential of the form 




(17) 



l/g{R). 




(18) 
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Fig. 4. Optimized liistograms for a dense two-dimensional Lennard- Jones fluid after 
feedback of the local diffusivity in a radial coordinate for varying temperatures. For 
the optimized ensemble peaks evolve at the free energy barriers between the shells 
of the liquid which proliferate at lower temperatures (solid lines). Additional peaks 
emerge at low temperatures (dashed lines) revealing interstitial states (arrows) 
between the shells of the liquid. 



where e is the interaction strength, a a length parameter, and i? the dis- 
tance between two particles. It is this distance i? between two arbitrarily 
chosen particles in the fluid that one can use as a new reaction coordinate 
for a projected random walk. For a given temperature defining an extended 
ensemble with weights Wcxtcndcd(^) and recording a histogram H(R) dur- 
ing a simulation will then allow to calculate the pair distribution function 
g{R) = i?(i?)/Wextendod(-R)- The pair distribution function g{R) is closely 
related to the potential of mean force (PMF) 

<?PMF(i?) - -^ln5(i?) , (19) 

which describes the average interaction between two particles in the fluid in 
the presence of many surrounding particles. 

For high particle densities and low enough temperatures shell structures 
form in the fluid which are reminiscent of the hexagonal lattice of the solid 
structure at very low temperatures. These shell structures are revealed by a si- 
nusoidal modulation in the PMF. Thermal equilibration between the shells is 
suppressed by entropic barriers which form between the shells. Again, one can 
ask what probability distribution, or histogram, should be sampled along the 
reaction coordinate, in this case the radial distance R, so that equilibration 
between the shells is improved. Measuring the local diffusivity for a random 
walk along the radial distance R in an interval [i?min, -f^max] and subsequently 
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applying the feedback algorithm described above optimized histograms H{R) 

arc found which arc plotted in Fig. 4 for varying temperatures [12]. The feed- 
back algorithm again shifts additional weight in the histogram towards the 
bottleneck of the simulation, in this case towards the barriers between the 
shells. Interestingly, additional peaks emerge in the optimized histogram as 
the temperature is lowered towards the vapor-liquid equilibrium. The minima 
between these peaks point to additional meta-stable configurations which oc- 
cur at these low temperatures, namely interstitial states which occur as the 
shells around two particles merge as detailed in Ref. [12]. 

This example illustrates that for some simulations the local diffusivity 
and optimized histogram them,selves are very sensitive measures that can 
reveal interesting underlying physical phenomena which are otherwise hard 
to detect in a numerical simulation. In general, a strong modulation of the 
local diffusivity for the random walk along a given reaction coordinate is 
a good indicator that the reaction coordinate itself is a good choice that 
captures some interesting physics of the problem. 

6 Pcirallel tempering / replica-exchange methods 

The simulation of frustrated and/or disordered systems suffers from a similar 

tunneling problem than the simulation of first-order phase transitions: local 
minima in energy space are separated by barriers that grow with system 
size. While the multicanonical or optimized ensembles do not help with the 
NP-hard problems faced by spin glasses, they are efficient in speeding up 
simulations of frustrated magnets without disorder [11]. 

An alternative to these extended ensembles for the simulation of frus- 
trated magnets is the "parallel tempering" or "replica-exchange" Monte Carlo 
method [13-16]. Instead of performing a single simulation at a fixed temper- 
ature, simulations are performed for M replicas at a set of temperatures 
Ti, T2, . . . , Tm- In addition to standard Monte Carlo updat(?s at a fixed tem- 
perature, exchange moves are proposed to swap two replicas between adjacent 
temperatures. These swaps are accepted with a probability 

min[l,exp(zi/3Z\£')], (20) 

where A(3 = (3j — Pi is the difference in inverse temperatures and AE 
I — Ei the difference in energy between the two replicas i and j. 

The effect of these exchange moves is that a replica can drift from a lo- 
cal free energy minimum at low temperatures to higher temperatures, where 
it is easier to cross energy barriers and equilibration is fast. Upon cooling 
(by another sequence of exchanges) it can end up in a different local mini- 
mum on time scales that are much shorter compared to a single simulation 
at a fixed low temperature. This random walk of a single replica in tem- 
perature space is the conjugate analog of the random walk in energy space 
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Fig. 5. Optimized temperature distribution H{T) for a two-dimensional Ising fer- 
romagnet with A'^ = 20 x 20 spins and 100 replicas/temperature points. After the 
feedback of the local diffusivity a peak evolves near the critical temperature of the 
transition Tc ~ 2.269. The feedback thereby shifts additional resources towards the 
bottleneck of the simulation which are identified by a suppressed local diffusivity. 
Note the similarity to Fig. 2. 

discussed for the extended ensemble techniques. The complement of the sta- 
tistical ensemble, defined by the weights Woxtondcd(£^), is the particular choice 
of temperature points in the temperature set {Ti, T2, . . . , Tm} for the parallel 
tempering simulation. The probability of sampling any given temperature T 
in an interval Ti < T < Tj+i can then be approximated by H[T) oc 1/ZiT, 
where AT — T^+i — Ti is the length of the temperature interval around the 
temperature T . This probability distribution H{T) is the equivalent to the 
histogram H{E) in the extended ensemble simulations. The ensemble op- 
timization technique discussed above can thus be reformulated to optimize 
the temperature set in a parallel-tempering simulation in such a way that 
the rate of round trips between the two extremal temperatures, Ti and Tm 
respectively, is maximized [17,18]. 

Starting with an initial temperature set {Ti, T2, . . . , Tm} a parallel tem- 
pering simulation is performed where each replica is labeled either "plus" or 
"minus" indicating which of the two extremal temperatures the respective 
replica has visited most recently. This allows to measure a current of replicas 
diffusing from the highest to the lowest temperature by recording two his- 
tograms, h^{T) and h^[T) for each temperature point. The current j is then 
given by 



J ^ D{T)H{T)-i- , 



(21) 



Ensemble optimization techniques 13 



4-x 



critical temperature 
T =2.269 

C 

X xjoodiiwocx X X 

I 
I 



M 
o 

•a 



xii^wii i:x X 



X XXX^IMKX XX X 



0- 



>no(xxx X X X ^ X 



4 6 

temperature T 



10 



Fig. 6. Optimized temperature sets for a two-dimensional Ising ferromagnet with 
= 20 X 20 spins. The initial temperature set with 20 temperature points is 
determined by a geometric progression [17,19] for the temperature interval [0.1, 10]. 
After feedback of the local diffusivity the temperature points accumulate near the 
critical temperature Tc = 2.269 of the phase transition (dashed line). Similar to the 
ensemble optimization in energy space the feedback of the local diffusivity relocates 
resources towards the bottleneck of the simulation. 



where D{T) is the local diffusivity for the random walk in temperature space, 
and f(T) = h+{T)/[h+(T) + h_(T)] is the fraction of random walkers which 
have visited the highest temperature Tm most recently. The probability dis- 
tribution H(T) is normalized, that is 



4 zt=^' 



where C is a normalization constant. Rearranging equation (21) the local 
diffusivity D{T) of the random walk in temperature space can be estimated 
as 

In analogy to the argument for the extended ensemble in energy space the 
current j is maximized by choosing a probability distribution 



iJopt.nal(r) OC - V ^ 1^ (24) 
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which is inversely proportional to the square root of the local diffusivity. The 

optimized temperature set {Tl,T2, . . . ,T^} is then found by choosing the 
n-th temperature point such that 



where M is the number of temperature points in the original temperature 
set, and the two extremal temperatures T{ = Ti and Tj^ = Tm remain 
unchanged. Similarly to the algorithm for the ensemble optimization this 
feedback of the local diffusivity should be iterated until the temperature set 
is converged. 

Figures 5 and 6 illustrate the so-optimized temperature sets for the Ising 
ferromagnet obtained by several iterations of the above feedback loop. After 
the feedback of the local diffusivity, temperature points accumulate near the 
critical temperature Tc = 2.269 of the transition. This is in full analogy to the 
optimized histograms for the extended ensemble simulations where resources 
are shifted towards the critical energy of the transition, for comparison see 
Figs. 2 and 5. 

It is interesting to note that for the so-optimized temperature set the 
acceptance rates for swap moves are not independent of the temperature 
[17]. Around the critical temperature, where temperature points are accu- 
mulated by the feedback algorithm, the acceptance rates are higher than at 
higher/lower temperatures, where the density of temperature points becomes 
considerably smaller after feedback [17,18]. The almost Markovian scaling be- 
havior for the optimized random walks in either energy or temperature space 
is thus generated by a problem-specific statistical ensemble which is char- 
acterized neither by a flat histogram nor flat acceptance rates for exchange 
moves, but by a characteristic probability distribution which concentrates 
resources at the minima of the measured local diffusivity. 



The ensemble optimization technique reviewed in this chapter should be 
broadly applicable to a wide range of applications possibly speeding up 
existing uniform sampling techniques by orders of magnitude. It has recently 
been used to improve broad-histogram Monte Carlo techniques [11] as well 
as parallel-tempering ^lontc Carlo simulations [17], with applications to frus- 
trated and disordered spin systems [11,17], dense fluids [12], as well as folded 
proteins [18]. It also holds promise to improve the simulation of quantum 
systems close to a phase transition when optimizing the {^xtcnded ensemble 
introduced for the quantum Wang-Landau algorithm outlined in Ref. [20]. 




(25) 



7 Outlook 
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